#include <iostream>
#include <opencv2/core.hpp>
#include <opencv2/imgcodecs.hpp>
#include <opencv2/imgproc.hpp>

static double _calc_correlation(const cv::Mat& hist1, const cv::Mat& hist2) {
    int n = hist1.rows * hist1.cols;
    double s1 = 0, s2 = 0;
    for (int i = 0; i < hist1.rows; i++) {
        for (int j = 0; j < hist1.cols; j++) {
            double f1 = hist1.ptr<float>(i)[j];
            double f2 = hist2.ptr<float>(i)[j];
            s1 += f1;
            s2 += f2;
        }
    }
    
    double m1 = s1 / n;
    double m2 = s2 / n;
    double s11 = 0, s12 = 0, s22 = 0;
    for (int i = 0; i < hist1.rows; i++) {
        for (int j = 0; j < hist1.cols; j++) {
            double f1 = hist1.ptr<float>(i)[j];
            double f2 = hist2.ptr<float>(i)[j];
            s11 += (f1-m1)*(f1-m1);
            s12 += (f2-m2)*(f1-m1);
            s22 += (f2-m2)*(f2-m2);
        }
    }
    double d = 0;
    double s11xs22 = s11 * s22;
    if (std::abs(s11xs22) > DBL_EPSILON) {
        d = s12 / std::sqrt(s11xs22);
    } else {
        d = 1.;
    }
    return d;
}

static double _calc_chi_square(const cv::Mat& hist1, const cv::Mat& hist2) {
    double d = 0;
    for (int i = 0; i < hist1.rows; i++) {
        for (int j = 0; j < hist1.cols; j++) {
            double f1 = hist1.ptr<float>(i)[j];
            double f2 = hist2.ptr<float>(i)[j];
            if (std::fabs(f1) > DBL_EPSILON) {
                d += (f1-f2)*(f1-f2)/f1;
            }
        }
    }
    return d;
}

static double _calc_intersection(const cv::Mat& hist1, const cv::Mat& hist2) {
    double d = 0;
    for (int i = 0; i < hist1.rows; i++) {
        for (int j = 0; j < hist1.cols; j++) {
            double f1 = hist1.ptr<float>(i)[j];
            double f2 = hist2.ptr<float>(i)[j];
            d += std::min(f1, f2);
        }
    }
    return d;
}

static double _calc_bhattacharyya_distance(const cv::Mat& hist1, const cv::Mat& hist2) {
    int n = hist1.rows * hist1.cols;
    double s1 = 0, s2 = 0, sqrt_s12 = 0;
    for (int i = 0; i < hist1.rows; i++) {
        for (int j = 0; j < hist1.cols; j++) {
            double f1 = hist1.ptr<float>(i)[j];
            double f2 = hist2.ptr<float>(i)[j];
            s1 += f1;
            s2 += f2;
            sqrt_s12 += std::sqrt(f1*f2);
        }
    }
    double s = s1*s2;
    s = (std::fabs(s) > FLT_EPSILON ? 1. / std::sqrt(s) : 1.);
    s = 1. - s*sqrt_s12;
    double d = std::sqrt(std::max(s, 0.));
    return d;
}

static double _calc_hellinger_distance(const cv::Mat& hist1, const cv::Mat& hist2) {
    return _calc_bhattacharyya_distance(hist1, hist2);
}

static double _calc_alternative_chi_square(const cv::Mat& hist1, const cv::Mat& hist2) {
    double d = 0;
    for (int i = 0; i < hist1.rows; i++) {
        for (int j = 0; j < hist1.cols; j++) {
            double f1 = hist1.ptr<float>(i)[j];
            double f2 = hist2.ptr<float>(i)[j];
            if (std::fabs(f1+f2) > DBL_EPSILON) {
                d += (f1-f2)*(f1-f2) / (f1+f2);
            }
        }
    }
    return d*2;
}

static double _calc_kullback_leibler_divergence(const cv::Mat& hist1, const cv::Mat& hist2) {
    double d = 0;
    for (int i = 0; i < hist1.rows; i++) {
        for (int j = 0; j < hist1.cols; j++) {
            double f1 = hist1.ptr<float>(i)[j];
            double f2 = hist2.ptr<float>(i)[j];
            if (std::fabs(f1) <= DBL_EPSILON) {
                continue;
            }
            if (std::fabs(f2) <= DBL_EPSILON) {
                f2 = 1e-10;
            }
            double dd = f1*std::log(f1/f2);
            d += dd;
        }
    }
    return d;
}

int main(int argc, char **argv) {
    cv::Mat src1 = cv::imread("/home/xlll/Downloads/opencv/samples/data/aloeL.jpg", cv::IMREAD_GRAYSCALE);
    cv::Mat src2 = cv::imread("/home/xlll/Downloads/opencv/samples/data/aloeR.jpg", cv::IMREAD_GRAYSCALE);
    if (src1.empty() || src2.empty()) {
        std::cout << "failed to read image" << std::endl;
        return EXIT_FAILURE;
    }
    
    int nimages = 1;
    int channels[] = {0};
    int dims = 1;
    int histsize[] = {256};
    float histrange[] = {0, 256};
    const float* ranges[] = {histrange};
    bool uniform = true;
    bool accumulate = false;
    
    cv::Mat hist1, hist2;
    cv::calcHist(&src1, nimages, channels, cv::Mat(), hist1, dims, histsize, ranges, uniform, accumulate);
    cv::calcHist(&src2, nimages, channels, cv::Mat(), hist2, dims, histsize, ranges, uniform, accumulate);
    
    double d1 = cv::compareHist(hist1, hist2, cv::HISTCMP_CORREL);
    double d11 = _calc_correlation(hist1, hist2);
    
    double d2 = cv::compareHist(hist1, hist2, cv::HISTCMP_CHISQR);
    double d22 = _calc_chi_square(hist1, hist2);
    
    double d3 = cv::compareHist(hist1, hist2, cv::HISTCMP_INTERSECT);
    double d33 = _calc_intersection(hist1, hist2);
    
    double d4 = cv::compareHist(hist1, hist2, cv::HISTCMP_BHATTACHARYYA);
    double d5 = cv::compareHist(hist1, hist2, cv::HISTCMP_HELLINGER);
    double d44 = _calc_bhattacharyya_distance(hist1, hist2);
    double d55 = _calc_hellinger_distance(hist1, hist2);
    
    double d6 = cv::compareHist(hist1, hist2, cv::HISTCMP_CHISQR_ALT);
    double d66 = _calc_alternative_chi_square(hist1, hist2);
    
    double d7 = cv::compareHist(hist1, hist2, cv::HISTCMP_KL_DIV);
    double d77 = _calc_kullback_leibler_divergence(hist1, hist2);
    
    std::cout << "d1 = " << d1 << std::endl;
    std::cout << "d11 = " << d11 << std::endl << std::endl;
    
    std::cout << "d2 = " << d2 << std::endl;
    std::cout << "d22 = " << d22 << std::endl << std::endl;
    
    std::cout << "d3 = " << d3 << std::endl;
    std::cout << "d33 = " << d33 << std::endl << std::endl;
    
    std::cout << "d4 = " << d4 << std::endl;
    std::cout << "d5 = " << d5 << std::endl;
    std::cout << "d44 = " << d44 << std::endl;
    std::cout << "d55 = " << d55 << std::endl << std::endl;
    
    std::cout << "d6 = " << d6 << std::endl;
    std::cout << "d66 = " << d66 << std::endl << std::endl;
    
    std::cout << "d7 = " << d7 << std::endl;
    std::cout << "d77 = " << d77 << std::endl;
    
    
    return 0;
}
